# ****************************************************************** #
# Part of the replication files for:

# Are Foreign Firms Good for the Environment? FDI and Protected Areas
# Ana Carolina Garriga and Muzhou Zhang

# See 0-README.txt for more information
# ****************************************************************** #



# R version 4.4.1 (2024-06-14)
# Last Updated Feb 14, 2025



# 0 Preparation -----

library(tidyverse)
library(cowplot)
library(haven)
library(sjlabelled)
library(psych)
library(xtable)

setwd("~/Replication-Garriga-Zhang-II") # CHANGE PATH TO WHERE THE FOLDER IS LOCATED IN YOUR MACHINE



# 1 Figure A1: Plot of Dependent Variable ----

data <- read_dta("1-data.dta") %>%
  mutate(country = sjlabelled::as_character(country))

ggplot(data) + geom_line(aes(year, papcnt)) +
  scale_x_continuous(limits = c(1980, 2020)) +
  facet_wrap(vars(country), ncol = 6, nrow = 10) +
  labs(x = "", y = "") +
  theme_classic() +
  theme(
    strip.text = element_text(size = 15),
    axis.text.x = element_text(angle = 45, size = 12, margin = margin(15, 0, 0, 0)),
    axis.text.y = element_text(size = 12),
    panel.grid.major.y = element_line(colour = "grey", linetype = "solid")
  )

ggsave("5-figures/supp-1-twoway.png", width = 15, height = 15)



# 2 Table A1: Summary Statistics ----

data_2summarize <- data %>%
  select(papcnt, fdistock, gdppc, kofgi, polity, popdens, cbd, diffusion, time)

describe(data_2summarize) %>%
  as_tibble(.) %>%
  mutate(Min = min,
         Median = median,
         Max = max,
         Mean = mean,
         SD = sd,
         Obs = n,
         .keep = "none") %>%
  add_column(Variable = get_label(data_2summarize), .before = 1) %>%
  xtable(caption = "Summary Statistics of Main Variables",
         align = c("l", "l", rep("r", 6)),
         digits = c(0, 0, rep(2, 5), 0))



# 3 Figure A2: LOCO Coefficients Plot ----

est_loco <- read_csv("3-loco-estimates.csv") %>% 
  slice(2:5) %>%
  mutate(X1 = c("b", "b_se", "a", "a_se")) %>%
  pivot_longer(loco1:loco60) %>%
  mutate(
    name,
    coef = if_else(str_detect(X1, "b"), "FDI stock", "FDI stock, squared (absolute value)"),
    stats = if_else(str_detect(X1, "se"), "se", "beta"),
    value = as.double(value),
    .keep = "none"
  ) %>%
  pivot_wider(
    names_from = "stats",
    values_from = "value",
    values_fn = list(value = abs)
  ) %>%
  mutate(
    across(beta:se, ~ as.double(.)),
    cilow = beta - 1.65 * se,
    cihigh = beta + 1.65 * se,
  )

ggplot(est_loco) +
  geom_pointrange(aes(y = name, x = beta, xmin = cilow, xmax = cihigh)) +
  facet_wrap(vars(coef), scales = "free", ncol = 2) +
  geom_vline(xintercept = 0, linetype = "longdash") +
  scale_x_continuous(expand = expansion(mult = c(0.1, 0.2))) +
  labs(x = "", y = "") +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    strip.background = element_rect(color = "black", fill = NA),
    panel.grid.major.x = element_line(colour = "grey", linetype = "dashed"),
    axis.line.y = element_blank(),
    axis.text.y.left = element_blank(),
    axis.ticks.y.left = element_blank(),
    aspect.ratio = 2/3,
    text = element_text(size = 20)
  )

ggsave("5-figures/supp-2-loco.png", width = 12, height = 8)



# **********************
# END OF R SCRIPT FILE # 
# **********************